Compare observed and model dune topography
Femke Batsleer
14 april 2021
#Load data
Field data, with predictions with INLA; interpolations with INLA
data from IBM model
#packages
library(tidyr)
library(dplyr)
library(INLA)
library(plotly)
library(forecast)
library(gridExtra)
#field data and predictions from INLA
lidar_data <- read.csv("output/Data_and_predictions_INLA_fielddata.csv", sep=",")
interlidar_data <- read.csv("output/interpolations_INLA_fielddata.csv", sep=",")
#data IBM output
ibm_data <- read.csv("data/data_model.csv", sep=",") %>% slice(-836) %>% #delete wrong (double) line
filter(error < 0.4) %>% #filter outlier from error
mutate(rep = as.factor(rep))#Transform variables
To get comparable values/units and scale them the same
#summary(lidar_data)
#summary(ibm_data)
#divide JC in two, than comparable to field data
ibm_data <- ibm_data %>% mutate(JC2 = JC/2) %>%
#standardize IBM data: P is from 0 to 1, so ok; Hstd scaled JC2
mutate(Pstd = P, Hstd = scale(JC2, center=mean(JC2), scale=sd(JC2)))#Models for IBM data
#first, get data for interpolations
interibm_data <- tibble(Pstd=seq(0,1, by=0.01), H=seq(20,70,by=0.5)) %>%
tidyr::expand(Pstd, H) %>%
mutate(Hstd=scale(H, center=mean(ibm_data$JC2), scale=sd(ibm_data$JC2))) %>% #scale by IBM data, to get similar values
mutate(volume = NA, error = NA)#add empty values for volume and error
#select needed columns from ibm_data
ibm_all <- ibm_data %>% select(Pstd, Hstd, volume, error) %>%
#bind with interpol data to feed it to INLA
bind_rows(select(interibm_data, -H))
#INLA
#for volume
f0_diff <- volume ~ Pstd + Hstd + Pstd:Hstd +
I(Pstd^2) + I(Hstd^2) +Hstd:I(Pstd^2) + Pstd:I(Hstd^2)
M0_diff <- inla(f0_diff,
#control.compute = list(dic = TRUE, waic = TRUE),
control.predictor = list(compute = TRUE),
family = "gaussian", data = ibm_all)
summary(M0_diff)##
## Call:
## c("inla(formula = f0_diff, family = \"gaussian\", data = ibm_all, ", " control.predictor = list(compute = TRUE))")
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 0.6524 10.7607 1.1690 12.5822
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode
## (Intercept) 55.8184 6.6912 42.7096 55.8082 68.9719 55.7884
## Pstd 84.3977 20.4905 44.0965 84.4218 124.5331 84.4720
## Hstd -31.6104 5.7107 -42.8082 -31.6158 -20.3933 -31.6261
## I(Pstd^2) -132.8997 19.4652 -171.0398 -132.9258 -94.6450 -132.9768
## I(Hstd^2) -4.3092 4.6528 -13.4432 -4.3102 4.8221 -4.3118
## Pstd:Hstd 68.0558 19.5353 29.6698 68.0658 106.3521 68.0877
## Hstd:I(Pstd^2) -12.2823 19.5982 -50.7301 -12.2932 26.1917 -12.3136
## Pstd:I(Hstd^2) 9.5700 7.8981 -5.9401 9.5700 25.0652 9.5708
## kld
## (Intercept) 0
## Pstd 0
## Hstd 0
## I(Pstd^2) 0
## I(Hstd^2) 0
## Pstd:Hstd 0
## Hstd:I(Pstd^2) 0
## Pstd:I(Hstd^2) 0
##
## The model has no random effects
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant
## Precision for the Gaussian observations 3e-04 0.00 2e-04 3e-04
## 0.975quant mode
## Precision for the Gaussian observations 3e-04 3e-04
##
## Expected number of effective parameters(std dev): 6.337(0.0272)
## Number of equivalent replicates : 156.06
##
## Marginal log-Likelihood: -5541.16
## Posterior marginals for linear predictor and fitted values computed#for error
f0_err <- error ~ Pstd + Hstd + Pstd:Hstd +
I(Pstd^2) + I(Hstd^2) +Hstd:I(Pstd^2) + Pstd:I(Hstd^2)
M0_err <- inla(f0_err,
control.predictor = list(compute = TRUE),
family = "gamma", data = ibm_all)
summary(M0_err)##
## Call:
## "inla(formula = f0_err, family = \"gamma\", data = ibm_all, control.predictor = list(compute = TRUE))"
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 0.5360 61.3286 1.4010 63.2656
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept) -4.5497 0.0569 -4.6608 -4.5499 -4.4375 -4.5504 0
## Pstd 0.9001 0.2278 0.4518 0.9005 1.3460 0.9012 0
## Hstd -0.1399 0.0511 -0.2401 -0.1399 -0.0395 -0.1400 0
## I(Pstd^2) -1.2839 0.2134 -1.7022 -1.2842 -0.8644 -1.2847 0
## I(Hstd^2) 0.0284 0.0304 -0.0311 0.0283 0.0881 0.0282 0
## Pstd:Hstd 0.9552 0.2226 0.5177 0.9553 1.3915 0.9556 0
## Hstd:I(Pstd^2) -0.6835 0.2213 -1.1178 -0.6837 -0.2489 -0.6839 0
## Pstd:I(Hstd^2) 0.0038 0.0520 -0.0983 0.0038 0.1057 0.0038 0
##
## The model has no random effects
##
## Model hyperparameters:
## mean sd 0.025quant
## Precision parameter for the Gamma observations 7.021 0.3052 6.435
## 0.5quant 0.975quant mode
## Precision parameter for the Gamma observations 7.016 7.634 7.007
##
## Expected number of effective parameters(std dev): 8.042(0.0019)
## Number of equivalent replicates : 122.98
##
## Marginal log-Likelihood: 3981.09
## Posterior marginals for linear predictor and fitted values computedsave(M0_diff, file="output/M0_diff.rda")
save(M0_err, file="output/M0_err.rda")
#extract interpolation data from IBM models
M0_diff_inter <- M0_diff$summary.fitted.values %>%
slice(990:n()) %>% #select predictions from dataset that are not input data
select(mean, sd, `0.025quant`, `0.975quant`) %>% #select neaded values
#rename to get similar values as in interlidar_data
rename(Int_diff_mean = mean, Int_diff_sd = sd, Int_diff_0.025quant=`0.025quant`, Int_diff_0.975quant = `0.975quant`)
M0_err_inter <- M0_err$summary.fitted.values %>%
slice(990:n()) %>% #select predictions from dataset that are not input data
select(mean, sd, `0.025quant`, `0.975quant`) %>% #select neaded values
#rename to get similar values as in interlidar_data
rename(Int_err_mean = mean, Int_err_sd = sd, Int_err_0.025quant=`0.025quant`, Int_err_0.975quant = `0.975quant`)
#bind these interpolations/predictions to interibm_data
interibm_data <- interibm_data %>% bind_cols(M0_diff_inter) %>% bind_cols(M0_err_inter)#Effect sizes
Visualisation of the effect sizes/signs of the covariates
#load INLA models from lidar data
load("output/M1_diff.rda")
load("output/M1_err.rda")
#lidar data get coefficients from models
post_beta1_diff <- M1_diff$summary.fixed[, c("mean", "0.025quant", "0.975quant")]
post_beta1_err <- M1_err$summary.fixed[, c("mean", "0.025quant", "0.975quant")]
post_beta1_diff <- post_beta1_diff %>% tibble::rownames_to_column(var="WhichVar") %>%
mutate(WhichVar = factor(WhichVar, levels=c("meandVstd", "Diststd", "Pstd:I(Hstd^2)", "Hstd:I(Pstd^2)",
"I(Hstd^2)","I(Pstd^2)", "Pstd:Hstd", "Hstd", "Pstd","Intercept" ))) %>%
mutate(WhichVar = recode(WhichVar, "meandVstd"="meandV", "Diststd"="Dist",
"Pstd:I(Hstd^2)"="P:JC²", "Hstd:I(Pstd^2)"="JC:P²",
"I(Hstd^2)"="JC²","I(Pstd^2)"="P²", "Pstd:Hstd"="P:JC", "Hstd"="JC",
"Pstd"="P","(Intercept)"="(Intercept)" ))
post_beta1_err <- post_beta1_err %>% tibble::rownames_to_column(var="WhichVar") %>%
mutate(WhichVar = factor(WhichVar, levels=c("meandVstd", "Diststd", "Pstd:I(Hstd^2)", "Hstd:I(Pstd^2)",
"I(Hstd^2)","I(Pstd^2)", "Pstd:Hstd", "Hstd", "Pstd","Intercept" ))) %>%
mutate(WhichVar = recode(WhichVar, "meandVstd"="meandV", "Diststd"="Dist",
"Pstd:I(Hstd^2)"="P:JC²", "Hstd:I(Pstd^2)"="JC:P²",
"I(Hstd^2)"="JC²","I(Pstd^2)"="P²", "Pstd:Hstd"="P:JC", "Hstd"="JC",
"Pstd"="P","(Intercept)"="(Intercept)" ))
#ibm data get coefficients from models
post_beta0_diff <- M0_diff$summary.fixed[, c("mean", "0.025quant", "0.975quant")]
post_beta0_err <- M0_err$summary.fixed[, c("mean", "0.025quant", "0.975quant")]
post_beta0_diff <- post_beta0_diff %>% tibble::rownames_to_column(var="WhichVar") %>%
mutate(WhichVar = factor(WhichVar, levels=c("Pstd:I(Hstd^2)", "Hstd:I(Pstd^2)",
"I(Hstd^2)","I(Pstd^2)", "Pstd:Hstd", "Hstd", "Pstd","(Intercept)" ))) %>%
mutate(WhichVar = recode(WhichVar, "Pstd:I(Hstd^2)"="P:JC²", "Hstd:I(Pstd^2)"="JC:P²",
"I(Hstd^2)"="JC²","I(Pstd^2)"="P²", "Pstd:Hstd"="P:JC", "Hstd"="JC",
"Pstd"="P","(Intercept)"="(Intercept)" ))
post_beta0_err <- post_beta0_err %>% tibble::rownames_to_column(var="WhichVar") %>%
mutate(WhichVar = factor(WhichVar, levels=c("Pstd:I(Hstd^2)", "Hstd:I(Pstd^2)",
"I(Hstd^2)","I(Pstd^2)", "Pstd:Hstd", "Hstd", "Pstd","(Intercept)" ))) %>%
mutate(WhichVar = recode(WhichVar, "Pstd:I(Hstd^2)"="P:JC²", "Hstd:I(Pstd^2)"="JC:P²",
"I(Hstd^2)"="JC²","I(Pstd^2)"="P²", "Pstd:Hstd"="P:JC", "Hstd"="JC",
"Pstd"="P","(Intercept)"="(Intercept)" ))
#lidar data
p_diff_lidar <- ggplot(post_beta1_diff) +
geom_pointrange() +
geom_hline(yintercept = 0, linetype = 2) +
aes(x=WhichVar, y =mean, ymin = `0.025quant`, ymax = `0.975quant`) +
xlab("covariate")+ ylab("Effect size")+
ggtitle("Posteriors for \u0394height (LiDAR data)")+
theme_classic() +
coord_flip()
#p_diff_lidar
p_err_lidar <- ggplot(post_beta1_err) +
geom_pointrange() +
geom_hline(yintercept = 0, linetype = 2) +
aes(x=WhichVar, y =mean, ymin = `0.025quant`, ymax = `0.975quant`) +
xlab("covariate")+ ylab("Effect size")+
ggtitle("Posteriors for CV(\u0394height) (LiDAR data)")+
theme_classic() +
coord_flip()
#p_err_lidar
#IBM model
p_diff_ibm <- ggplot(post_beta0_diff) +
geom_pointrange() +
geom_hline(yintercept = 0, linetype = 2) +
aes(x=WhichVar, y =mean, ymin = `0.025quant`, ymax = `0.975quant`) +
xlab("covariate")+ ylab("Effect size")+
ggtitle("Posteriors for volume (IBM data)")+
theme_classic() +
coord_flip()
#p_diff_ibm
p_err_ibm <- ggplot(post_beta0_err) +
geom_pointrange() +
geom_hline(yintercept = 0, linetype = 2) +
aes(x=WhichVar, y =mean, ymin = `0.025quant`, ymax = `0.975quant`) +
xlab("covariate")+ ylab("Effect size")+
ggtitle("Posteriors for CV(volume) (IBM data)")+
theme_classic() +
coord_flip()
#p_err_ibm
grid.arrange(p_diff_ibm, p_err_ibm, p_diff_lidar, p_err_lidar)
#Visualisation of the interpolations
3D graphs for both the model and lidar data
Difference (height/volume): lijkt me bij combinate hoge P en H en kleine P en H enkel verschillend van elkaar. Bij lidar data gaat bij hoge P en H gaat difference ferm naar beneden, bij IBM data niet. Bij lage P en H gaat bij ibm data de difference naar omhoog, bij lidar data weer naar beneden. Bij de niet extreme waarden lijkt het wel goed overeen te komen (zou eens apart kunnen geplot worden).
Error: zadelprofiel bij ibm data, bij lidar data niet aanwezig (grote en kleine waarden van H)…
#3Dplot
#Lidar data
Interlidar_diff_3D <- plot_ly(data = interlidar_data, x = ~H, y = ~Pstd, z=~Int_diff_mean, color=~Int_diff_mean) %>% add_markers() %>%
layout(title = "Difference in dune height (LiDAR data)",
scene = list(
xaxis = list(title = "JC"),
yaxis = list(title = "P", tickvals=c(0,0.2,0.4,0.6,0.8,1)),
zaxis = list(title = "\u0394height (m)", titleangle=90))) %>%
colorbar(title = "Prediction \u0394height")## Warning: `arrange_()` is deprecated as of dplyr 0.7.0.
## Please use `arrange()` instead.
## See vignette('programming') for more help
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.Interlidar_err_3D <- plot_ly(data = interlidar_data, x = ~H, y = ~Pstd, z=~exp(Int_err_mean), color=~exp(Int_err_mean)) %>% add_markers() %>%
layout(title = "Variability in dune height (LiDAR data)",
scene = list(
xaxis = list(title = "JC"),
yaxis = list(title = "P", tickvals=c(0,0.2,0.4,0.6,0.8,1)),
zaxis = list(title = "CV (m)"))) %>%
colorbar(title = "Prediction CV")
#IBM model
Interibm_diff_3D <- plot_ly(data = interibm_data, x = ~H, y = ~Pstd, z=~Int_diff_mean, color=~Int_diff_mean) %>% add_markers() %>%
layout(title = "Difference in dune volume (IBM data)",
scene = list(
xaxis = list(title = "JC"),
yaxis = list(title = "P", tickvals=c(0,0.2,0.4,0.6,0.8,1)),
zaxis = list(title = "volume (m³)"))) %>%
colorbar(title = "Prediction volume")
Interibm_err_3D <- plot_ly(data = interibm_data, x = ~H, y = ~Pstd, z=~exp(Int_err_mean), color=~exp(Int_err_mean)) %>% add_markers() %>%
layout(title = "Variability in dune volume (IBM data)",
scene = list(
xaxis = list(title = "JC"),
yaxis = list(title = "P", tickvals=c(0,0.2,0.4,0.6,0.8,1)),
zaxis = list(title = "CV (m³)"))) %>%
colorbar(title = "Prediction CV")
Interlidar_diff_3DInteribm_diff_3DInterlidar_err_3DInteribm_err_3D#Plotting uncertainties
#3Dplot
#Lidar data
Interlidar_diff_3D_uncertainty <- plot_ly(data = interlidar_data, x = ~H, y = ~Pstd, z=~Int_diff_mean, color=~(Int_diff_0.975quant - Int_diff_0.025quant)) %>% add_markers() %>%
layout(title = "Uncertainty Difference (lidar data)",
scene = list(
xaxis = list(title = "JC"),
yaxis = list(title = "P", tickvals=c(0,0.2,0.4,0.6,0.8,1)),
zaxis = list(title = "Height (m)")))
Interlidar_err_3D_uncertainty <- plot_ly(data = interlidar_data, x = ~H, y = ~Pstd, z=~exp(Int_err_mean), color=~exp(Int_err_0.975quant-Int_err_0.025quant)) %>% add_markers() %>%
layout(title = "Uncertainty Variability (lidar data)",
scene = list(
xaxis = list(title = "JC"),
yaxis = list(title = "P"),
zaxis = list(title = "Error (m)")))
Interibm_diff_3D_uncertainty <- plot_ly(data = interibm_data, x = ~H, y = ~Pstd, z=~Int_diff_mean, color=~(Int_diff_0.975quant - Int_diff_0.025quant)) %>% add_markers() %>%
layout(title = "Uncertainty volume (ibm data)",
scene = list(
xaxis = list(title = "JC"),
yaxis = list(title = "P"),
zaxis = list(title = "Height (m)")))
Interibm_err_3D_uncertainty <- plot_ly(data = interibm_data, x = ~H, y = ~Pstd, z=~exp(Int_err_mean), color=~exp(Int_err_0.975quant-Int_err_0.025quant)) %>% add_markers() %>%
layout(title = "Uncertainty variability (ibm data)",
scene = list(
xaxis = list(title = "JC"),
yaxis = list(title = "P"),
zaxis = list(title = "Error (m)")))
Interlidar_diff_3D_uncertaintyInteribm_diff_3D_uncertaintyInterlidar_err_3D_uncertaintyInteribm_err_3D_uncertainty